﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "aitken")]
    public static unsafe double aitken(IntPtr x_ptr, IntPtr y_ptr, int n, double eps, double t)
    {
        double* x = (double*)x_ptr.ToPointer();
        double* y = (double*)y_ptr.ToPointer();

        return aitken(x, y, n, eps, t);
    }

    /// <summary>
    /// Aitken逐步插值
    /// </summary>
    /// <param name="x">x[n]存放n个给定的有序结点值。</param>
    /// <param name="y">y[n]存放n个给定结点上的函数值。</param>
    /// <param name="n">n给定结点的个数。</param>
    /// <param name="eps">eps插值的精度要求。</param>
    /// <param name="t">t指定插值点值。</param>
    /// <returns>函数返回指定插值点t处的函数近似值。</returns>
    public static unsafe double aitken(double* x, double* y, int n, double eps, double t)
    {
        int i, j, k, m, l = 0;
        double z;
        double[] xx = new double[10];
        double[] yy = new double[10];

        z = 0.0;
        if (n < 1) return (z);
        if (n == 1)
        {
            z = y[0];
            return (z);
        }
        //最多取前后10个点
        m = 10;
        if (m >= n) m = n;
        //起始点
        if (t <= x[0]) k = 1;
        //起始点
        else if (t >= x[n - 1]) k = n;
        else
        {
            k = 1; j = n;
            while ((k - j != 1) && (k - j != -1))
            {
                l = (k + j) / 2;
                if (t < x[l - 1]) j = l;
                else k = l;
            }
            //起始点
            if (Math.Abs(t - x[l - 1]) > Math.Abs(t - x[j - 1]))
            {
                k = j;
            }
        }
        j = 1;
        l = 0;
        //从起始点开始轮流在前后取m个点
        for (i = 1; i <= m; i++)
        {
            k = k + j * l;
            if ((k < 1) || (k > n))
            {
                l = l + 1;
                j = -j;
                k = k + j * l;
            }
            xx[i - 1] = x[k - 1];
            yy[i - 1] = y[k - 1];
            l = l + 1;
            j = -j;
        }
        i = 0;
        //对m个点作Aitken逐步插值
        do
        {
            i = i + 1; z = yy[i];
            for (j = 0; j <= i - 1; j++)
            {
                z = yy[j] + (t - xx[j]) * (yy[j] - z) / (xx[j] - xx[i]);
            }
            yy[i] = z;
        } while ((i != m - 1) && (Math.Abs(yy[i] - yy[i - 1]) > eps));
        return (z);
    }

    /*
    // Aitken逐步插值例
      int main()
      { 
          double t,z,eps;
          double x[11]={-1.0,-0.8,-0.65,-0.4,-0.3,
                             0.0,0.2,0.4,0.6,0.8,1.0};
          double y[11]={0.0384615,0.0588236,0.0864865,0.2,
                0.307692,1.0,0.5,0.2,0.1,0.0588236,0.0384615};
          eps=1.0e-6;
          t=-0.75; z=aitken(x,y,11,eps,t);
          cout <<"t = " <<t <<"     z = " <<z <<endl;
          t=0.05; z=aitken(x,y,11,eps,t);
          cout <<"t = " <<t <<"     z = " <<z <<endl;
          return 0;
      }
    */
}

